UpdateSoilMoisture Subroutine

SUBROUTINE UpdateSoilMoisture  & 
!
(id, i, j, dt, soilDepthCell, percolationcellRZ, percolationcellTZ, &
 runoffcell, Qin, Qout )

IMPLICIT NONE

! Arguments with intent(in):
INTEGER (KIND = short), INTENT(IN) :: id !!soil balance id
INTEGER (KIND = short), INTENT(IN) :: i,j
INTEGER (KIND = short), INTENT(IN) :: dt
REAL (KIND = double), INTENT(IN) :: soilDepthCell
REAL (KIND = float), INTENT(IN) :: Qin
REAL (KIND = float), INTENT(IN) :: Qout



!Arguments with intent(inout):
REAL (KIND = double), INTENT(INOUT) :: runoffcell
REAL (KIND = float), INTENT(INOUT) :: percolationcellRZ
REAL (KIND = float), INTENT(INOUT) :: percolationcellTZ

!Local declaration:
REAL (KIND = double) :: smPrevRZ !! root zone soil moisture of the previous time step
REAL (KIND = double) :: smPrevTZ !! transmission zone soil moisture of the previous time step
REAL (KIND = float)  :: check
REAL (KIND = double) :: dsRZ !!root zone soil mositure excess rate (m/s)
REAL (KIND = double) :: dsTZ !!transmission zone soil mositure excess rate (m/s)
REAL (KIND = float)  :: sdRZ  !!root zone soil depth (m)
REAL (KIND = float)  :: sdTZ  !!transmission zone soil depth (m)
REAL (KIND = double) :: subIn, subOut  !!subsurface fluxes
REAL (KIND = double) :: area !!cell area (m2)
!------------end of declaration------------------------------------------------

IF ( id == LANDPLAIN .AND. saturatedByGroundwater ) THEN
    ! water table depth lies within the root zone or above ground
    percolationcellRZ = 0.
    percolationcellTZ = 0.
    infilt % mat(i,j) = 0.
    et % mat(i,j) = pet % mat(i,j)
    capRise % mat (i,j) = et % mat(i,j)
    runoffcell = rainBalance % mat (i, j)
    soilMoistureRZ % mat(i,j) = thetas % mat (i,j)
    soilMoistureTZ % mat(i,j) = thetas % mat (i,j)
    soilMoisture % mat(i,j) = thetas % mat (i,j)
    balanceError % mat (i,j) = 0.
    RETURN
END IF

IF ( id == LAKE ) THEN
    soilMoistureRZ % mat(i,j) = thetas % mat (i,j)
    soilMoistureTZ % mat(i,j) = thetas % mat (i,j)
    soilMoisture % mat(i,j) = thetas % mat (i,j)
    balanceError % mat (i,j) = 0.
    deltaSoilMoisture % mat (i,j) = 0.
    RETURN
END IF

smPrevRZ = soilMoistureRZ % mat(i,j)
smPrevTZ = soilMoistureTZ % mat(i,j)

sdRZ = soilDepthRZ % mat (i,j)
sdTZ = soilDepthTZ % mat (i,j)

!IF ( soilDepthCell < soilDepthRZ % mat(i,j) ) THEN
!    sdRZ = soilDepthCell
!    sdTZ = 0.
!ELSE
!    sdRZ = soilDepthRZ % mat(i,j)
!    sdTZ = MAX (0.1, ( soilDepthCell - sdRZ ) ) 
!END IF

IF ( id == HILLSLOPE .OR. id == CHANNEL ) THEN
    area = CellArea (soilMoisture, i, j) 
    subIn = Qin / area
    subOut = Qout / area
ELSE IF ( id == LANDPLAIN ) THEN
    subIn = 0.
    subOut = 0.
END IF


IF (sdTZ > 0.) THEN
	soilMoistureTZ % mat(i,j) = smPrevTZ + & 
                        ( percolationcellRZ -  &
                        percolationcellTZ + &
                        subIn - subOut      )  /  &
                        sdTZ * dt
        
    IF (soilMoistureTZ % mat(i,j) < thetar % mat(i,j)) THEN
        soilMoistureTZ % mat(i,j) = thetar % mat(i,j)
    END IF
        
    IF ( soilMoistureTZ % mat(i,j) >  thetas % mat(i,j) ) THEN !saturation excess
        dsTZ = ( soilMoistureTZ % mat(i,j) - thetas % mat(i,j) ) * sdTZ / dt
            soilMoistureTZ % mat(i,j) = thetas % mat(i,j)
    ELSE
        dsTZ = 0.
    END IF
              
END IF

IF (sdRZ > 0.) THEN
	soilMoistureRZ % mat(i,j) = smPrevRZ + ( infilt % mat(i,j) - percolationcellRZ - &
                            et % mat(i,j) + capRise % mat(i,j) + dsTZ  )  / sdRZ * dt		   
              
END IF


!particular cases
!case 1: soil is too dry, soil moisture goes below residual value
IF (soilMoistureRZ % mat(i,j) < thetar % mat(i,j) ) THEN
  !compute soil moisture deficit (negative) converted to flux (m/s)
  dsRZ = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) * soilDepthCell / dt
  
  !adjust evapotranspiration
  et % mat(i,j) = et % mat(i,j) + dsRZ
     
  !check for negative et
  IF (et % mat(i,j) < 0. ) et % mat(i,j) = 0.
     
  !limit soil moisture
  soilMoistureRZ % mat(i,j)  = thetar % mat(i,j)
  
  
!case 2: soil is too wet, soil moisture goes above saturated value  
ELSE IF ( soilMoistureRZ % mat(i,j) > thetas % mat(i,j) ) THEN
  !compute soil moisture excess converted to flux (m/s)
	dsRZ = ( soilMoistureRZ % mat(i,j) - thetas % mat(i,j) ) * sdRZ / dt
	!adjust runoff for saturation excess runoff and infiltration rate
	runoffcell = runoffcell + dsRZ
	infilt % mat(i,j) = infilt % mat(i,j) - dsRZ
    
    !limit soil moisture to saturation value
    soilMoistureRZ % mat(i,j)  = thetas % mat(i,j)
  
    !check if infiltration is negative after correction
    !adjust capillary rise if present (plain cell)
	IF (infilt % mat(i,j) < 0.) THEN
      dsRZ = - infilt % mat(i,j)
      infilt % mat(i,j) = 0.
      IF (id == LANDPLAIN ) THEN
         capRise % mat(i,j) = capRise % mat(i,j) - dsRZ
      END IF
	END IF
	

END IF

!balance error
check = infilt % mat(i,j) +   &
        capRise % mat(i,j) -  &
        percolationcellTZ  -    &
        et % mat(i,j) +       &
        subIn - subOut - &
        (soilMoistureRZ % mat(i,j) - smPrevRZ) * sdRZ / dt - &
        (soilMoistureTZ % mat(i,j) - smPrevTZ) * sdTZ / dt

balanceError % mat (i,j) = check * dt * 1000.


!soil mositure variation
deltaSoilMoisture % mat (i,j) = (soilMoistureRZ % mat(i,j) - smPrevRZ) * sdRZ + &
        (soilMoistureTZ % mat(i,j) - smPrevTZ) * sdTZ 


!update relative saturation and mean soil moisture
soilSatRZ % mat(i,j) = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) / &
                     ( thetas % mat(i,j) - thetar % mat(i,j) )
soilSatTZ % mat(i,j) = ( soilMoistureTZ % mat(i,j) - thetar % mat(i,j) ) / &
                     ( thetas % mat(i,j) - thetar % mat(i,j) )
soilSat % mat(i,j) = ( soilSatRZ % mat(i,j) * sdRZ  + &
                       soilSatTZ % mat(i,j) * sdTZ  ) / (sdRZ + sdTZ)
IF ( sdRZ + sdTZ > 0. ) THEN
     soilMoisture % mat(i,j) = ( soilMoistureRZ % mat(i,j) * sdRZ + &
                            soilMoistureTZ % mat(i,j) * sdTZ ) / (sdRZ + sdTZ)
ELSE
    soilMoisture % mat(i,j) = thetas % mat (i,j) 
END IF

RETURN
END SUBROUTINE UpdateSoilMoisture